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CO Abstract 

We prove that any stable method for resolving the Gibbs phenomenon - that is, recovering 
high-order accuracy from the first m Fourier coefficients of an analytic and nonperiodic function 
- can converge at best root-exponentially fast in m. Any method with faster convergence must 
also be unstable, and in particular, exponential convergence implies exponential ill-conditioning. 
This result is analogous to a recent theorem of Platte, Trefethen & Kuijlaars concerning recovery 
from pointwise function values on an equispaced m-grid. The main step in our proof is an 
estimate for the maocimal behaviour of a polynomial of degree n with bounded m-term Fourier 
series. This bound is related to a conjecture of Hrycak & Grochenig. In the second part of the 
paper we discuss the implications of our main theorem to polynomial-based interpolation and 
^ least-squares approaches for overcoming the Gibbs phenomenon. Finally, we propose the use of 

so-called Fourier extensions as an attractive alternative for this problem. We present numerical 
results demonstrating rapid convergence of the resulting approximation in a stable manner. 
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1 Introduction 

The Fourier series 




of an analytic and periodic function / : [— 1, 1] — > M converges exponentially fast in the truncation 
parameter m. However, such rapid convergence is destroyed once periodicity is no longer present. 
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For nonperiodic /, the series J-mf{x) converges only linearly in compact subsets of (—1, 1), and there 
is no uniform convergence on [—1, 1]. Near the endpoints x = ±1 one witnesses the well-known Gibbs 
phenomenon [34j[49]. 

Although the Gibbs phenomenon has a long history [35], and is completely understood mathe- 
matically |341 149] ■ it remains a significant hurdle in practical applications of Fourier series. Indeed, 
it is a testament to its importance that the question of the resolution of the Gibbs phenomenon - i.e. 
restoring high-order convergence given only the partial Fourier series J>„/ of a function - remains 
an area of active inquiry. 

Although there are many existing methods for this problem (see Section [2] for a review), no one 
stands out as being inherently superior. In particular, many exponentially convergent strategies 
appear to suffer from some sort of ill-conditioning, and whilst there are other methods that resolve 
the Gibbs phenomenon in a stable manner, these typically result in slower orders of convergence. 
The purpose of this paper is to explain these observations. 

Specifically, we show that any exponentially-convergent method must also be exponentially ill- 
conditioned, and in general, if a method has a convergence rate of p~™ for some p > 1 and 
T e (^,1], then it must also possess ill-conditioning of order p™ . In particular, this result 
implies the following fundamental stability barrier: the best possible convergence rate for a stable 
method for resolving the Gibbs phenomenon is root-exponential in m. 

A theorem of this type is not new. Our main result is a direct analogue of a theorem of Platte, 
Trefethen & Kuijlaars for the problem of overcoming the Runge phenomenon in equispaced poly- 
nomial interpolation. In |45j it was proved that any method for recovering high accuracy from the 
pointwise values of a function on an equispaced grid must also exhibit the aforementioned instability 
behaviour. The problem we consider in this paper, namely recovering high accuracy from the first 
2m + 1 Fourier coefficients of /, can be considered a continuous analogue of this problem. Indeed, 
the problem of recovery from equispaced pointwise values is equivalent to that of recovery from 
discrete Fourier coefficients. 

In proving our main result we follow a similar argument to that of [45^. The key step therein 
is the use of an estimate of Coppersmith & Rivlin [21] concerning the maximal behaviour of an 
algebraic polynomial of degree n which is bounded on an equispaced grid of m points. Our main 
theoretical contribution is an analogous result for Fourier series: namely, we estimate the maximal 
behaviour of a polynomial p of degree n with bounded m-term Fourier series J-mP- In doing so, we 
provide a partial answer to a conjecture of Hrycak & Grochenig [321 (see Section [s] for a discussion). 

Our main theorem on the behaviour of polynomials with bounded Fourier sums has an impor- 
tant consequence for a particular method that is sometimes used in practice. Our result implies 
that the so-called inverse polynomial reconstruction method (IPRM) |4H 1421 144j is exponentially 
unstable. Moreover, although it is possible to stabilize this method via a least squares procedure 
(henceforth referred to as the polynomial least squares method), our main theorem demonstrates 
that this necessarily decreases the convergence rate to root-exponential. 

Although root-exponential convergence is the best possible permitted, our theorem says nothing 
about spectral convergence. Nor does is apply to methods for which convergence occurs only down 
to some finite tolerance etoi- In the final part of this paper we propose the use of so-called Fourier 
extensions [3 [HI [HI HB HO] for this problem. As we demonstrate by numerical example, this approach 
gives rapid convergence - sometimes exponential, but always spectral - but only down to a finite 
tolerance on the order of machine precision. Moreover, we show that it typically outperforms the 
aforementioned polynomial least squares method. We conclude that Fourier extensions present an 
attractive method for this problem. 

The outline of the remainder of this paper is as follows. In Section [2] we give an overview 
of common methods for the resolution of the Gibbs phenomenon, and discuss the key issues of 
convergence and stability. The aforementioned estimate is obtained in Section [3| and in Section [4] 
we present the main result of the paper. In Section [sj we consider polynomial interpolation and least 
squares, and Fourier extensions are introduced in Section [6] 
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2 The Gibbs phenomenon and its resolution 



The Gibbs phenomenon has a long history dating back to Wilbraham in 1848 [SD] (to acknowledge 
the contribution of Wilbraham, the name Gibbs- Wilbraham phenomenon is also occasionally used). 
Forgotten for half a century, this phenomenon was rediscovered by Michelson |43j , with the ensuing 
debate regarding convergence, or lack thereof, between Michelson and Love (carried out in Nature) 
being eventually settled by Gibbs [351 ES] in 1899, after the arbitration of Poincare. Gibbs contri- 
bution to this problem was first recognised by Bocher in 1906 |13j . who introduced the term the 
Gibbs phenomenon. A detailed and fascinating review of the Gibbs phenomenon and its history is 
provided in [38J , with shorter summaries also presented in \22\ I34j . 

Many methods have been proposed to ameliorate or resolve the Gibbs phenomenon. Of these, 
perhaps the earliest to appear were filters and moUifiers. Here the Gibbs phenomenon is viewed as 
noise polluting the high-order Fourier coefficients, which can therefore be mitigated by premultipli- 
cation with a rapidly decaying function [311 HH] • Unfortunately, standard filters do not lead to high 
uniform accuracy: they only ensure faster convergence in regions of [—1, 1] away from the endpoints 
X = ±1. More recently, attention has focused on adaptive filters and moUifiers (see [49] and refer- 
ences therein), which can be constructed to obtain exponential convergence in (—1,1) in a stable 
manner, with typically polynomial accuracy at the endpoints. Note that this does not contradict 
the main result of this paper, since the rate of uniform convergence is not exponential. As a general 
principle, exponential convergence away from x = ±1 can be obtained without ill-conditioning. 

An alternative approach (which can also be viewed as a type of moUifier [49]) is the technique of 
spectral reprojection, introduced and developed by Gottlieb et al [351 [331 |331|351|3S] . Here the slowly 
convergent Fourier series is reprojected onto a suitable basis; the so-called Gibbs complementary basis. 
Often Gegenbauer polynomials are used for this basis. Provided the various parameters are carefully 
tuned, exponential convergence, uniformly on [—1,1], can be restored. Unfortunately, the original 
Gegenbauer reconstruction procedure of [35] has been shown to be rather sensitive to the choice 
of parameters [171 I28j . with the wrong parameters giving potentially divergent approximations. 
To mitigate this effect, a substantially more robust procedure, based on Freud polynomials, was 
introduced in [25]. Nonetheless, our main result states that spectral reprojection must either exhibit 
exponential instability or not be truly exponentially convergent. We note at this point that a stability 
analysis of spectral reprojection has not yet been carried out. Nevertheless, this approach has found 
application in a number of areas, including image processing [9] llOj and the spectral approximation 
of PDEs with discontinuous solutions [37] . 

Spectral reprojection is sometimes referred to a direct technique, since it does not require solution 
of a linear system. The most obvious inverse method is to seek to 'interpolate' the first 2m -I- 1 
Fourier coefficients of / with an algebraic polynomial of degree 2m + 1. This technique, which 
requires solution of a linear system of equations, is sometimes referred to as the inverse polynomial 
reconstruction method (IPRM) in literature |411l42lHl] . However, interpolation of Fourier coefficients 
can be seen as a continuous analogue of polynomial interpolation on equispaced nodes. It should come 
as little surprise, therefore, that there are substantial issues with both convergence and stability. In 
particular, a Runge-type phenomenon is witnessed. See ^[55], as well as Sectionjsj for a discussion. 

Since 'interpolating' Fourier coefficients may not work, one can also use a lower degree polynomial 
in combination with a least squares fit (so-called polynomial least squares). This was first discussed 
in detail in ^39j, and later in \4^. Unfortunately, as we shall later prove, to ensure stability one 
requires that the degree n scale like -^/m. This corresponds to only root-exponential convergence in 
m, consistently with the stability barrier we establish in Section [4] Nonetheless, despite this slower 
convergence, we remark that this approach often outperforms spectral reprojection in practice. For 
a comparison, see [4]. 

As an alternative to lowering the polynomial degree, one may also try using a higher degree 
polynomial. Underdetermined least squares leads to a poor approximation in this case. However, 
better accuracy can be restored by using /^-minimization instead. To the best of knowledge no 
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analysis currently exists for this approach. For a related discussion in the case of equispaced function 
values, see |19LI45) . One may also consider Sobolev norm minimization, such has been considered in 
the equispaced case in [23] , 

In 0] a general framework was introduced for stable reconstructions in Hilbert spaces. Given 
Fourier coefRcients, one can reconstruct in any other basis of functions, with one example being the 
polynomial least squares method discussed above. An alternative to polynomials involves the use 
of splines. As discussed in [S] , fixed-order splines result in algebraic convergence in a stable manner 
(see also (51]). One may also consider variable-order splines, but stability becomes an issue. 

A different approach to overcoming the Gibbs phenomenon is to smooth the function / via 
subtraction so as to make it periodic up to a given order, and then compute its Fourier series. This 
idea dates back to Krylov and Lanczos, amongst others, and was later studied by Lax and Gottlieb 
& Orszag - see [TJ HI HB] and references therein. Such smoothing can be carried out implicitly, 
via extrapolation on the high-order Fourier coefficients; an approach sometimes known as Eckhoff 's 
method in literature [26]. As shown in [46], this method converges algebraically fast in to. However, 
there are also issues with instability [1]. A hybrid approach, combining Gegenbauer reconstruction 
and Eckhoff 's method, was also developed in [27] . 

Alternative methods for overcoming the Gibbs phenomenon arise from sequence extrapolation 
techniques. See [12 [2D]. In a similar spirit, DriscoU & Fornberg introduced Fade-based method in 
[25j . This approach gives exponential convergence [llj . however in view of our theorem, must also 
be exponentially unstable. Such instability was noted in [25j . 

There are numerous other methods for resolving the Gibbs phenomenon, and we have not pre- 
sented a complete list. The reader is referred to [THl [IH] and references therein for further 
information. We also mention that many methods designed for the related problem of recovering 
high accuracy from function values on equispaced grids (i.e. overcoming the Runge phenomenon) 
can potentially be adapted to the Fourier coefficient problem. See [19l [45] for a comprehensive list 
of such methods. 

One such technique that has recently been successfully applied to overcome the Runge phe- 
nomenon is that of Fourier extensions. In the final part of this paper (Section [6]) we propose the use 
of this approach to overcome the Gibbs phenomenon. As we show, the corresponding method can 
be extremely effective for this problem. 

3 Maximal behaviour of an algebraic polynomial with bounded 
Fourier series 

The first step towards the stability theorem in [45] is an estimate due to Coppersmith & Rivlin [21] . 
This concerns the behaviour of the quantity 



An,rr. = SUp : P G P„, = l} , 



where ||p||^ = sup^g[_i i] \p{x)\, \\p\\^^^ = maxj=i,._„ |p(a;j)| and {xi,...,Xm} is a grid of m 
equispaced points in [—1,1]. In [21] it was shown that 

(ci) " < ^„,™ < (cs)"^. (3.1) 
for constants C2 > ci > 1. This estimate determines how many equispaced gridpoints are required 



to control the behaviour of a polynomial of degree n. Observe that if m = o(n^) then (3.1 1 implies 
that there exists a polynomial which is bounded on the grid {xi, . . . , Xm\^ but which grows large in 
between grid points. On the other hand, An^m is bounded as 7i,to — > oo if and only if to = O (n^) 
(that the scaling m ~ O (n^) is sufficient for boundcdness is an older result which dates back to 
Schonhage [48]). 
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The problem of the behaviour of An^m is a classical one in approximation theory (see [JS] for 
a summary of its history). However, A„.m involves equispaced grid values. To prove a stability 
theorem for overcoming the Gibbs phenomenon, we need to study a related quantity involving 
Fourier coefficients. To this end, we define 

Bn,m = sup{|b||2 : p e P„, b||„, = 1} , (3.2) 



where ||p||2 = y j\ b(a;)p dx is the L^-norm on [—1, 1] and 



Ml - E \Pj\' = \\^mP\\\ (3.3) 

ljl<™ 

is the Z^-norm of the first 2m + 1 Fourier coefficients of p. 

We remark that the quantity i?„,m differs from An^m in two respects. First, it involves continuous 
Fourier coefficients, as opposed to discrete Fourier coefficients (i.e. equispaced grid values). Second, 
rather than using the uniform norm, we consider the (respectively P)-norm. This is natural, 
since Fourier series satisfy Parseval's relation in the L^/P-norms. Indeed, the second equality in 
(3.3) follows directly from this relation. 

With this aside, note that 

hm Bn,rn = 1, 

m— J-oo 

for any fixed n € N. This follows from strong convergence of the operators J-m — > I on L^( — 1, 1), 
and the fact that the space P„ is finite dimensional. Moreover, much as in the case of Aji^jji^ it can 
be shown that the scaling m — O (n^) is sufficient for boundedness of Bn,m- This result was first 
proved by Hrycak & Grochenig [3^1 Thm. 4.1] (see also [U Lem. 3.1]). 

Our main result in this section shows that Bn,m admits a similar lower bound to that found in 
(Isll). We have 



Theorem 3.1. Let B„^m be as in (3.2). Then there exists a constant c > 1 such that 

Bn,m > c^, \/n,m G N. 

Specifically, we have the lower bound 

{Bn,mf > 1 + ^ + TT^rf " , (3.4) 

where d > |. 

Proof. Let p € P„ be arbitrary. Integrating by parts, we find that 



where 



Therefore, we can write 



1 



V2(i7r)'= 

n 

= (-1)^^(7), 3 e Z\{0}, m E • (3.5) 



fc=i 
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Note that p £ P„ is uniquely defined by the values po, 61, ... , &„. Therefore, (3.5 1 defines a one-to- 
one correspondence between polynomials p = p{x) with pQ = 0, say, and polynomials p(t) satisfying 
p(0) = 0. Hence, using Parseval's relation, we have 

Bn,m > sup{||p||2 : p G Vn,Po = 0, = 1} = sup B{n,m,p), 

per* 

= {p e P„ : p{0) = 0} and 



where 



B{n, m,p) 



Si<|j|<oo b( j)P 

\ Z]l<|j|<m 



To prove the theorem it is sufficient to find a particular P e P* such that B{n,m, P) admits the 
bound Q. 

Let n = Aq + 1 for some g € N. Consider the polynomial 

P{x) ■.^xT*{x^)A,ix^), 

where 

Mx'):=f[ix'-],), 

and T* (x) denotes the Chebyshev polynomial of degree q on the interval 
definition, 



m2 ' (q+l)2 



Then, by 



and therefore we have 



(B{n,m,P)f 



l<|j|<oo 



\P{] 



1\\2 



g<|j"|<oo 



\P{] 



1\\2 



El 



'1^12 



^ 1 < I J?' I < m 

Note that Aq{-) has all its zeros in the interval [^,1], hence \Aq\ is monotonically decreasing on the 
interval [0, Therefore, 

< \Aq{j,)\ < \Aq{^)\ < , q<ji<m<j2. 



1^12 



1 



E 



7n<j<oo 



|P(} 



1M2 



E« 



(3.6) 



So, putting expression for P into (3.6), we obtain 



{B{n,m,P)f =1 + 



>1 



E 



=1 



Eq<j<m j^\Pq i pWl^qi pW 
^-^m<j'<C30 j2 

■iT;(^)nA,(;i,)p 

Eq<j<m J2|r*(j2-)p|^q(;;^)p 

V J_|7n*f'lM2 



rM2 



E<;<j<m j2|r*(^.2)h 



:=1 



Since \Tq\ < 1 on the interval [— 



m2 ' (q+l)' 



N 
D 

], for the denominator D we have the estimate 



D 



q<j<oo 
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so that D < i. As for the numerator N, we spht it into two parts: 



m<j<2m ' 



2m<j<oo 



By definition, T* is the Chebyshev polynomial of degree q on the interval ^^^1)2 ]; so all of 
its zeros are located inside that interval. Hence |T*| is decreasing on [0,^] towards the value 
= 1- Therefore 



Ni > 
N2 > 



m<j <2m 

^9 ((2^)1 E > 2k ^9 ((2^) 



2m<j<oo ' 



I.e. 



^ 2m ' ^ 2m 



1 



9 ^(2m)^ 



Let us evaluate the value |T*(a;5)|, where a;Q := j2rnj^' Take the afhne mapping 

M:{r = [^,j^]}^{I=[-l,l]}, 

so that 

T*{x) = T,{M{x)) ^ T;(4) = T,(M(x5)) = T,{x^) , 

where Tq{x) is the standard Chebyshev polynomial on / = [— 1, 1]. The length of the interval /* is 
less than ^ , so mapping M onto the interval / with the length 2 uses the length magnifying factor 

A > 2g^. The point Xq = j2m)^ ^^^^ distance = from the left endpoint ^ of /*, so it 

will be mapped to the point xq < —1 given by 

-a;o = l + <5o, 8o>2q^5*o = ^. 

For a; = 1 + (5, we have 

Tq{x) = \({x+ ^/x^ - If + {x- \/x2 - > \{x+^x'^-lY> \{l + V2~5Y, 
and since \/25o = \/3^ > we have 



9 /»" 



■ 2 ' 



Hence 



m 
9 ' 



(3.7) 
(3.8) 



Combing these results together, and using the fact that g « , we obtain 

N 



{Bin,m,P)r >l + ^ = l + ^ + § 

>1 + JL + JLj^lVm 

2m Am 

>1 + ^ + :^7"^/«™ • 
8m 16m 



Since m > n/2 and < n/4, we have r = ^ > 2 so that a lower bound for 7 is 7 = (1 + > 9/4. 



This completes the proof. 



□ 



7 



10" 
10'" 

10' 

1000- 



10' 

100- 



10' - 
10' 
1000 
100 

10- 



0.6- 
05 - 
0.4: 
03 - 



1.2 - 
1.0 - 

0.8; 

0.6 

0.4 - 
0.2 - 



(«,/?) = (ii) 



= (it) 



= (if) 



Figure 1: Top row: the quantity B. 
row: the scaled quantities n^~^ log(S, 
ematica using additional precision. Note that the quantity Bn,m can be computed, since concides with the 
minimum singular value of a particular matrix. 



n an/5 (squarcs) and the lower bound B* (circles) against n. Bottom 
n an?) n^~^ log(B* ^^fi) ■ Computatious were carried out in Math- 



In Figure [T] we confirm this theorem by plotting the quantity Bn,m and the lower bound 



b: 



n /9 



8m 16m V 4 



Note that B* ^ not only provides a lower bound, it also appears to correctly predict the behaviour 
of Bn^m- We conjecture that an upper bound of the form 

Bn^m<c'^, Vn,meN, (3.9) 
also holds, much as in the case of the quantity An^m- 

4 An impossibility theorem for the resolution of the Gibbs 
phenomenon 

We now turn our attention to the main result of this paper. First, we require some notation. 
Following [35], let {4>m}meN be a family of mapping 1, 1) — >• 1, 1) such that (t)m{f) depends 
only on the values {fj}\j\<m- Note that the mappings (j)m can be both hnear or nonlinear. We define 
the condition number for 0^ by 

Umif + 9) - cj^minh 

Km — sup hm sup 77— r. . 

/ 9; ^ llsIL 

0<llslL<e 

For a compact set £^ C C we shall also let B{E) be the Banach space of functions continuous on E 
and analytic in its interior, with norm ||/||^ = sup^^^ 
Our main result is as follows: 
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Theorem 4.1. Let a compact set E C C containing [—1, 1] in its interior be fixed, and suppose that 
{0m}mgN such that (i) for any f , the quantity (j^mif) depends only on the values {fj}\j\<mj f^i^d 
(a) for some M < oo, a > 1 and r g (|, I]? "we have 

\\f-4>,n{f)h<Ma~"'^\\f\\E, yfeB{E),meN. (4.1) 
Then the condition numbers satisfy 

for some c > 1 and all sufficiently large m. 

Proof. We proceed as in [45, . We may without loss of generahty assume that £^ is a Bernstein eUipse 



||/-0™(/)||2 < m>mo, yfeB{E{p)). (4.2) 



E{p) for some p> I. We may then replace ( |4.l[ ) by 

1 

2' 

where toq is sufhciently large and a > is fixed. Let p e P„ . An inequality of Bernstein |451 Lemma 
1] implies that 

\\p\\e<p1\p\\[-i,i]- 

Now, since p € P„, we have that < cri||p||2 for some constant c > 0. Thus, 

IblU <C7ip"|b||2. (4.3) 



Setting p = / in (4.2 1 now gives 

Um{p)h > IHI2 - Ip—^blU > (1 - icnp™^) \\ph. 
Suppose that n < ^ara^ . Then 

cnp"-""'" < icam^/9~5"™" < 1^ Vm > toi, 
where mi is sufficiently large. Set m* :~ max{r7io, mi}. If m > m* this now gives 

||</'m(p)||2 > 5lb||2, ym>m*. 
Since 0m (0) — (this follows from (|4.1[)) we have 



</>«(£?)-(?!>« (0)11 2 _ ||0m(ep)||2 ^ 1 ||ej3||2 _ 1 b|l2 



hP\L imL ~ 2 ||ep||„, 2 IIpII 

Therefore, since p G P,i is arbitrary, we obtain 



Km> h sup 114^ r = hB{n,m) 



Vn < |am^. 



where B{n,m) is given by (3.2 1. Using Theorem 3.1 with n ~ ^ani^ , this now yields 



> 1 ia^m^--^ 



as required. □ 

This theorem implies that any method for overcoming the Gibbs phenomenon that results in a 
convergence rate of cr^™ for all analytic functions / G B{E) must also possess ill-conditioning of 
order c™^^ \ In particular, the best convergence rate that can be achieved with a stable method 
is root-exponential in the number of Fourier coefficients m. As commented in [45] . we stress that, 
despite the use of polynomials in the proof, this theorem is not about polynomials or polynomial- 
based approximation procedures. It holds for all (linear or nonlinear) mappings ipm satisfying a 



bound of the form (4.1). 
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5 Fourier coefficient interpolation and least squares 



As discussed in Section [2] an obvious way to seek to overcome the Gibbs phenomenon is by inter- 
polating the first 2m + 1 Fourier coefhcients of / with a polynomial of degree 2m. In other words, 
we construct an approximation /,„ satisfying 



fmj = fj, |j| < m, /„ 



■ 2m- 



(5.1) 



It can be shown that such a polynomial exists uniquely for any m [39j. However, as we shall see 
in a moment, this approach is both exponentially unstable and divergent. A simple modification 
involves computing an overdetermined least squares with (n < to): 



n,m = argmm 




(5.2) 



Note that fm, as defined by (5.1), coincides with fn.m when n 



As mentioned, (5.11 is often 



referred to as the inverse polynomial reconstruction method (IPRM) [HI The modification 

(5.2), referred to as polynomial least squares, was introduced and analysed in [39], and developed 



further in [J] (see also [3]). 
In [B] it was shown that /„ 



I satisfies the sharp bound 

\\f - fn,m\\ < B2n,m inf 



where Bn,m is the constant defined in (3.2 1 (previous, but non-sharp, estimates were given in [H 139)1. 
Moreover, it was also shown that the condition number k ~ Hn.m of the mapping / i— > fn.m is 
precisely 



where Bn.m is defined by (3.2). Hence, Theorem 3.1 allows us to explain both the convergence and 

(5.4) 



stability of this approach. In particular, (5.3) and Theorem 3.1 give that 

^n,m ^ C ^ . 

and therefore polynomial least squares is unstable whenever n grows faster than y/rri. In the partic- 
ular case n = m, this shows that the IPRM method (5.1 ) is exponentially ill-conditioned. Note that 



such exponential growth was previously observed, although not analysed, in [51 1421 

We remark that Hrycak & Grochenig have conjectured that an upper bound of the form (5.4) 

71^ 

also holds. In other words, the condition numbers can grow no worse than for some c > 1. This 



is of course equivalent to the conjecture (3.9 1 concerning the constant Bn.m- 

Regarding the convergence of fn,m, it is useful to recall that the error infpgp^n 11/ ~ P\\ decays 
like where p > 1 is the parameter of the largest Bernstein ellipse within which / is analytic 

[37]. Thus, we have the bound 

11/ — fn,m\\ < CfBn.mP^^" , 



(5.5) 

where Cf > is some constant depending on / only. This indicates that fn,m may fail to converge to 
/ if the rate of growth of Bn^m exceeds that of p^". In particular, if n = 0{m), in which case Bn.m 
is exponentially large, there will always be functions (analytic within only a small ellipse E(p)) for 
which the right-hand side of (5.5) diverges. Thus, the polynomial least squares method ( |5.2[ ), and 
in particular the IPRM (5.1), may well suffer from a Runge-type phenomenon - i.e. divergence of 
fn.m for some nontrivial family of analytic functions - whenever n = O (to). Although we have no 



proof of this fact (the upper bound (5.5) need not be sharp for an individual function /), there is 
substantial numerical evidence to support this conjecture I39| 1421 . 

On the other hand, when n ^ O (a/to) it is known that Bn.m = O (1) [U[3n|, and therefore we have 
stability, as well as guaranteed convergence of the approximation. Unfortunately, the convergence 
rate is only root-exponential, and Theorem |4. 1| demonstrates that it can be no faster. 
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Remark 5.1 Overdetermined least squares is a well-known approach to overcome the Runge phe- 
nomenon in equispaced polynomial interpolation [141 145| . In |14) essentially the same arguments 
as those given above were presented for this problem, leading to the same conclusions: namely, a 
Runge- type phenomenon for n = 0{ni), but stability and root-exponential convergence whenever 
n = 0{y/m). 

We note that the use of polynomial least squares for overcoming the Gibbs phenomenon, as op- 
posed to the Runge phenomenon, is far less well known. However, this approach (with n — O {y/rn)) 
does appear to outperform other more commonly used techniques, such as spectral reprojection, 
despite its formally slower convergence [4j. 

6 Fourier extensions for overcoming the Gibbs phenomenon 

The principle behind the IPRM and polynomial least squares is to reconstruct a function / in a 
finite-dimensional subspace in which it is well-approximated, with the space of polynomials of degree 
n being a natural choice for analytic /. Recently these ideas were developed substantially in a series 
of paper by Adcock & Hansen [H [51 |B] . Therein a framework, known as generalized sampling, was 
introduced to stably reconstruct elements of abstract Hilbert spaces in finite-dimensional subspaces 
from their samples taken with respect to a particular basis or frame. Polynomial least squares is 
one specific instance of this general framework, based on sampling with respect to the Fourier basis 
and reconstructing in the subspace P2„. 

However, since generalized sampling allows reconstructions in arbitrary spaces, there is no need 
to choose this particular space. In this final section of this paper, we consider the use of an alter- 
native subspace for reconstruction, and demonstrate that this gives substantial improvements over 
polynomial least squares. 

The particular subspace we shall employ is based on an approximation scheme known as Fourier 
extensions. Fourier extensions have been used in the past to successfully overcome the Runge 
phenomenon. In particular, Boyd [161 119] and Bruno et al [21] both show that this approximation 
can recover analytic functions to extremely high accuracy from equispaced data. This was confirmed 
by the analysis of Huybrechs [40], and later Adcock et al [8], wherein it was shown that Fourier 
extensions based on equispaced function data converge exponentially fast for all functions analytic 
in a sufficiently large region, but only down to an error on the order of machine precision. Hence there 
is no contradiction with the impossibility theorem of [45| (which requires exponential convergence 
for all m). For functions analytic within only small regions, the convergence is spectral down to 
approximately the same level. 

The idea behind Fourier extensions is very simple. We seek to approximate an analytic and 
nonperiodic function / : [— 1, 1] M using a Fourier series on an extended domain [— T, T], where 
T > 1 is fixed. In other words, we compute an approximation from the finite-dimensional subspace 



Note the construction fn,m is almost identical to that considered in [5J dni IHI] for the case of 
equispaced data, the only difference being that the discrete least-squares is taken over a sum of 
Fourier coefficients as opposed to pointwise function values. Note also that fn,m is similar to the 




For the problem considered in this paper, where only the Fourier samples {fj}\j\<m of / are given, 
we propose the following Fourier extension approximation to /: 




(6.1) 
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polynomial least squares approximation (5.2). However, critically the approximation space has been 
changed from the polynomial space P2n to the Fourier extension space iS„. 

For the sake of brevity, we shall not analyse fn,m- Instead, we show by numerical example how 
effective this approximation can be in practice. We remark in passing that the analysis of can be 
adapted the case of fn,m,- In particular, one can verify identical convergence behaviour of /„.,„ as 
in the case of equispaced data. 



To demonstrate the effectiveness of (6.1) we compare it to the polynomial least squares ap- 
proximation /„,,„ given by (5.2). Recall that both methods involve additional parameters: n, the 
polynomial degree, in the case of /„,m, and n, the degree of the Fourier extension, and T, the size of 
the extension interval, in the case of fn,m- Hence, given to, we need some way to select these values 
in order to make a comparison. We shall do this as follows. In the case of fn,m we let n be as large 
as possible whilst keeping the condition number kpls < '*0: where kq is some fixed value (we use 
kq = 10 in our experiments), and for fn.m we first fix T and then proceed in the same way so that 
'*FE l£ I'^o- In other words, for both approximations we ensure that the condition number is no worse 
than kq. Thus both methods are guaranteed to be equally robust with respect to perturbations. 

Note that computing the kpls is straightforward: provided an orthonormal polynomial basis 
of scaled Legendre polynomials is used, one only needs to determine the minimal singular value 
of a particular matrix [4^. On the other hand, computing the condition number for the Fourier 
extension is more challenging. As was explained in [8j, one cannot determine the condition number 



by simply examining singular values of the matrix of the linear system resulting from (6.1 ). Instead 



as proposed in we compute the condition number as follows. Since the method is linear, we have 

\\Fn,rrM\2 

KpE = sup r-j , 

fc#0 

where Fnrn{b) := ^^S^'^'^{^\^\<m\^3 ~ '/'jP} is the Fourier extension computed from the vector 

of Fourier coefficients h. We can therefore approximate kfe by randomly selecting vectors h and 
computing ||F„.,„(6)||2/||6||;2. Hence, if b^^\ . . . S C^™"*"^ are chosen uniformly at random with 
||6[^l||/2 < 1, j = 1, . . . , we consider the approximation 

\\Fn.n{b^'^)h 

KPE = . max — f^T- w KFE. 

In all our experiments we take the number of trials t — 100. In order to compute Kpg we also need 
to approximate ||i^„,m(^'''')||2- We do this with an equispaced quadrature based on 2001 nodes. 

In Figure [2] we plot for each method the computed values of n against the number of samples 
TO. Note that n = O {y/rn) for polynomial least squares, exactly as the results of the last section 
predict. On the other hand n grows linearly in to for the Fourier extension, with the constant of 
proportionality depending on the choice of the parameter T. 

Using these values, in Figure [3] we plot the errors committed by the two methods for a variety 
of different test functions. We consider four types of functions: 

(i) Entire functions with boundary layers, e.g. f{x) — e"^^^^), a ^ 1. 

(ii) Meromorphic functions with real singularities near a; = 1, e.g. f{x) — , a S> 1. 

(iii) Meromorphic functions with complex singularities on the imaginary axis near x = 0, e.g. the 
classical Runge example f{x) = ^^^2^2 , a ^ 1. 

(iv) Entire, oscillatory functions, e.g. f{x) = coswTra;, w ^ 1. 
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0.2 



50 100 150 200 50 100 150 200 

polynomial least squares Fourier extension 

Figure 2: Top row: The parameter n against m = 1, . . . ,200, where n is chosen as large as possible so 
that < 10 for the Fourier extension and kpls < 10 for polynomial least squares. In the right column, 
squares, circles and crosses correspond to the parameter values T = |,2,4 respectively. Bottom row: the 
scaled parameter n/^/rn (left) and n/m (right). 



Note that algebraic polynomials are particularly well suited to (i) and (ii). One can resolve a 
boundary layer of width 1/a using a polynomial of degree O {^/a) [151 131j . Fourier approximations 
on the other hand, such as the Fourier extension, are less well suited for boundary layers, since they 
require O (a) degrees of freedom. However, with only Fourier data prescribed, the polynomial least 
squares method can stably reconstruct a polynomial of degree at most 0{yjm), meaning that a 
boundary layer of width a cannot be resolved any more efficiently in terms of the number of Fourier 
coefficients m by polynomial least squares than by the Fourier extension. Hence, this particular 
advantage of using a polynomial reconstruction space disappears. As can be seen in Figure [3| fn,m 
and fn,m give roughly the same errors for functions of type (i) and (ii). 

On the other hand, Fourier extensions significantly outperform polynomial least squares for 
functions of type (iii) and (iv). This can be explained in a similar manner. Up to constant factors, 
Fourier extensions and algebraic polynomial approximations require the same number of degrees of 
freedom n to resolve functions of type (iii) and (iv) [7]. However, since we recover a Fourier extension 
of degree O (m), whereas we can only recover a polynomial of O {^/rn), we obtain a vastly superior 
approximation using the former. 

Notice also one interesting facet of Figure [3j namely, the choice of T has little effect on the ap- 
proximation error. This is due to the balance between degrees of freedom (i.e. n) and approximation 
properties. Larger T means a larger n can be used for a given m whilst retaining the condition num- 
ber (see Figure [2|. However, larger T also translates into slower (although still spectral/exponential) 
convergence in n 0- Perhaps surprisingly, these two effects balance in practice, rendering the choice 
of T insignificant. To the best of our knowledge, this observation has not previously been made in 
the literature on Fourier extensions. 

In summary, Fourier extensions appear to present an extremely attractive method for the problem 
of recovering analytic functions from Fourier data. A more thorough comparison, incorporating more 
of the methods discussed in Section [2| is the topic of future work. 
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Figure 3: Errors for the polynomial least squares and Fourier extension approximations fn,m and fn,m 
against m — 10, 20, . . . , 200. Squares correspond to polynomial least squares, circles, crosses and diamonds 
correspond to Fourier extensions with parameter T = §,2,4 respectively. For each m, the parameter n was 
determined so that the condition number of the corresponding method was at most 10. 
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